function [sigma_draw, beta_draw] = BVARdrawpost(VAR)
% =======================================================================
% Draw from the distribution of beta and sigma
% =======================================================================
% [sigma_draw, beta_draw] = BVARdrawpost(VAR)
% -----------------------------------------------------------------------
% INPUT
%   VAR  : structure, result of VARmodel function
% OUPUT
%   sigma_draw : draw from the posterior of VCV matrix of residuals of VAR
%   beta_draw  : draw from the posterior of beta of VAR
% =======================================================================
% Ambrogio Cesa Bianchi, July 2011
% ambrogio.cesabianchi@gmail.com



%% Get relevant parameters from VAR structure
%============================================
nobs = VAR.nobs;
k     = [];
nvars = VAR.neqs;
X     = VAR.X;

% likelihood estimates
beta_hat    = VAR.beta;
sigma_hat   = VAR.sigma;
inv_sig_hat = inv(sigma_hat);
    
%% Draw the VCV matrix (sigma)
%=============================
inv_sigma_draw = wishrnd(inv_sig_hat/nobs, nobs);
sigma_draw     = inv(inv_sigma_draw);

%% Draw the coeffiecient matrix (beta)
%=====================================
con_3     = inv(X'*X);
con_4     = kron(sigma_draw, con_3);
Bhat_vec  = beta_hat(:);
draw      = mvnrnd(Bhat_vec,con_4);
beta_draw = reshape(draw, k, nvars);
